Breakfast food preferences of 42 subjects for 15 breakfast foods. The data are the rank orders of preference (1 = most prefered). Here are the first five rows and 9 columns of the data set.
library("smacof")
library("rgl")
breakfast[1:5, 1:9]
## toast butoast engmuff jdonut cintoast bluemuff hrolls toastmarm butoastj
## 1 13 12 7 3 5 4 8 11 10
## 2 15 11 6 3 10 5 14 8 9
## 3 15 10 12 14 3 2 9 8 7
## 4 6 14 11 3 7 8 12 10 9
## 5 15 9 6 14 13 2 12 8 7
Compute 2D and 3D unfolding solutions.
res2D <- smacofRect(breakfast, ndim = 2)
res3D <- smacofRect(breakfast, ndim = 3)
res4D <- smacofRect(breakfast, ndim = 4, itmax = 5000)
res2D
##
## Call: smacofRect(delta = breakfast, ndim = 2)
##
## Model: Rectangular smacof
## Number of subjects: 42
## Number of objects: 15
##
## Stress-1 value: 0.324
## Number of iterations: 114
res3D
##
## Call: smacofRect(delta = breakfast, ndim = 3)
##
## Model: Rectangular smacof
## Number of subjects: 42
## Number of objects: 15
##
## Stress-1 value: 0.296
## Number of iterations: 546
res4D
##
## Call: smacofRect(delta = breakfast, ndim = 4, itmax = 5000)
##
## Model: Rectangular smacof
## Number of subjects: 42
## Number of objects: 15
##
## Stress-1 value: 0.282
## Number of iterations: 1329
op <- par(pty = "s")
mxplt <- max(abs(max(res2D$conf.row, res2D$conf.col)), abs(min(res2D$conf.row, res2D$conf.col))) * 1.1
plot(res2D, type = "p",
pch = 19,
xlim = c(-mxplt, mxplt),
ylim = c(-mxplt, mxplt),
col.columns = "red",
col.rows = "lightblue",
label.conf.columns = list(label = TRUE, pos = 3, col = "red"),
label.conf.rows = list(label = TRUE, pos = 3, col = "lightblue"))
abline(v = 0, h = 0, col = "gray")
par(op)
mxplt <- max(abs(max(res3D$conf.row, res3D$conf.col)), abs(min(res3D$conf.row, res3D$conf.col))) * 1.1
plot3d(res3D$conf.col, col = "red",
type = "s",
radius = .3,
xlim = c(-mxplt, mxplt),
ylim = c(-mxplt, mxplt),
zlim = c(-mxplt, mxplt))
text3d(res3D$conf.col,
text = rownames(res3D$conf.col),
cex = 0.7,
adj = c(1.6, 1.6))
# add the subjects' ideal points
plot3d(res3D$conf.row, col = "lightblue",
add = TRUE,
type = "s",
radius = 0.2)
text3d(res3D$conf.row,
text = rownames(res3D$conf.row),
cex = 0.5,
adj = c(1.5, 1.5))